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SUMMARY 



three-dimensional time dependent Navier- Stokes analysis was applied to the 
rotor blade vortex interaction problem. The numerical procedure is an iterative 
implicit procedure using three point central differences to represent spatial 
derivatives. A series of calculations were made to determine the time steps, 
pseudo- time steps, iterations, artificial dissipation level, etc. required to 
maintain a nondissipative vortex. Results show the chosen method to have 
excellent non -dissipative properties provided the correct parameters are chosen. 
This study was used to set parameters for both two- and three-dimensional blade 
vortex interaction studies. The two-dimensional study considered the interaction 
between a vortex and a NACA0012 airfoil. The results showed the detailed physics 
during the interaction including the pressure pulse propagating from the blade. 
The simulated flow physics was qualitatively similar to that experimentally 
observed. The 2-D BVI phenomena is the result of the buildup and violent 
collapse of the shock waves and local supersonic pockets on the blade surfaces. 
The resulting pressure pulse build-up appears to be centered at the blade leading 
edge. The three-dimensional interaction study considered the case of a vortex at 
20° incidence to the blade leading edge. Although the qualitative results were 
similar to that of the two-dimensional interaction, details clearly showed the 
three-dimensional nature of the interaction process. 


INTRODUCTION 


The interaction ofv concentrated vortices with blades induces unsteady 


aerodynamic loading responsible for blade vibrations, aeroelastic instabilities, 


% 




.V 


and impulsive noise. The effects of blade-vortex interaction (BVI) are especially 

significant in the transonic flow K regime, in which the strength and position of 

■ \ 

the shock waves are sensitive to small changes in the flow parameters. Due to its 

commonv occurrence in many aerodynamic applications, BVI has been a subject of many 

N 

experimental, analytic and computational investigations. Recent examples of 

\ "\ 

experimental investigations include the work\of Caradonna and his colleagues 

\ 

(1984, 1988) and Booth and Yu (1986). Analytical studies include those of Panares 
(1987), Lee and Smith (1987) and Poling, Dadone and Telionis (1989). 

Investigations based upon numerical simulations include the two-dimensional full 
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Navier-Stokes simulations of Sankar and Tang (1985), the two-dimensional thin 
layer Navier-Stokes simulations of Srinivasan and McCroskey (1987), Srinivasan, 
McCroskey and Baeder (1986) and Rai (1987), and the three-dimensional thin layer 
Navier-Stokes simulation of Srinivasan and McCroskey (1989). 

The BVI problem can be viewed as an unsteady, three-dimensional close 
encounter of curved vortex filaments, at arbitrary intersection angles, with a 
blade that is in combined translational and rotational motion. Under certain 
operating conditions, the blade can encounter a vortex that is almost parallel to 
itself. Such an encounter is essentially two-dimensional but unsteady, and has 
been the focus of many BVI investigations. However, in general the blade vortex 
interaction is a three-dimensional phenomenon. Both two-dimensional and three- 
dimensional interactions are considered here. 

At the present time, a key problem in computing flows containing concentrated 
vortices is the ability to preserve and convect these vortices in a finite - 
difference or finite-volume grid without false numerical diffusion due to 
truncation error, artificial dissipation and turbulence modelling. Various 
investigators have dealt with the vortex preservation problem in different ways. 
One approach, the prescribed vortex approach (e.g. Srinivasan and McCroskey, 1987 
and 1989) assumes a background flow field consisting of a specified vortex and 
solves the governing partial differential equations for the difference between the 
full dependent variable field and that of the specified vortex. Under this 
approach the numerical diffusion is applied only to that part of the flow field 
which represents the perturbation from the isolated vortex flow field. A 
conceptually similar approach is taken by Sankar and Tang (1985) who solve the 
Navier-Stokes equations for the usual dependent variables; i.e., do not use a 
perturbation approach but modify the numerical dissipative operator to act on the 
difference between the instantaneous total field values and some presumed vortex 
field values, thereby removing a large part of the spurious dissipation of the 
vortex structure. More recently, strong BVI problems were solved by using a 
fifth-order accurate upwind-biased scheme without employing any vortex 
preservation techniques (Rai, 1987). 

In the present work, the ensemble -averaged, time -dependent Navier-Stokes 
equations are solved on a body- fitted grid around a NACA0012 airfoil in two and 
three dimensions to study strong interaction of a vortex with a stationary blade. 
The Navier-Stokes equations are solved by. using an iterative implicit finite- 
difference scheme with second order spatial and temporal accuracies. Furthermore, 
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simple vortex preservation techniques are used to minimize the amount of spurious 
numerical dissipation and eddy viscosity caused by the presence of the vortex 
during its convective motion towards the leading edge of the blade. In terms of 
spatial differencing, the present study employs three -point central differences 
with adjustable dissipation terms. An alternative approach is to use various 
upwind differencing schemes in transformed coordinates. For example, Rai (1987) 
has used first- and second- order upwind differences and also fifth- order 
upwind-biased differences to study blade vortex interactions . Although 
higher -order centered and upwind schemes can be more accurate, the order of 
accuracy refers to transformed coordinate variables. The high order schemes are 
usually first-order accurate in physical space variables for general non-uniform 
grids. Centered schemes can be made equivalent to upwind schemes by using special 
forms of dissipation with particular choices of dissipation parameters, as 
discussed in Pulliam (1985), Yoon and Kwak (1988) and Briley, Govindan and 
McDonald (1990) . 

One attraction of the upwind schems is that they propagate waves according to 
local characteristic behavior (in a one -dimensional sense). This is offset 
somewhat by complexity and the use of expanded grid molecules for higher order 
upwind schemes, which complicates their use near boundaries. An upwind approach 
can also be used to advantage in constructing special schemes for capturing strong 
shocks. Central-difference schemes also approximate local physics correctly for 
smooth flows, but require special treatment of (extraneous) boundary conditions to 
avoid error due to boundary conditions which conflict with characteristic behavior 
at boundaries . 

In our experience, the present scheme is more accurate than a first-order 
upwind scheme, and has accuracy comparable to the higher-order upwind schemes for 
general non-uniform grids. The present scheme also has a three -point bandwidth 
and is simple and economical. The simple artificial dissipation approach used 
here does not require assumptions regarding the velocity field associated with the 
prescribed vortex approach used by some other investigators, e.g. Srinivasan et 
al. (1989). Finally, the main focus of the present effort is a three-dimensional 
simulation based upon the full Navier- Stokes equations whereas the previous 
three-dimensional simulations, Srinivasan et al. (1989), are based upon a solution 
of the thin layer equations. 

The present effort proceeds in three separate steps. First consideration is 
given to the isolated vortex to develop a run protocol which has suitable vortex 
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preservation properties. Vortex preservation is necessary for viable two- and/or 
three-dimensional simulations. Secondly, a two-dimensional, transonic blade 
vortex interaction simulation is made. This allows a qualitative assessment prior 
to the three-dimensional interaction. The third portion of the program is the 
three-dimensional blade vortex interaction simulation . 


II. ANALYSIS 

Governing Equations and Solution Algorithm 

Prior to discussing the results a discussion of the governing equations and 
numerical approach is appropriate. As previously stated, the present approach 
solves the full ensemble -averaged Navier-Stokes equations without any shear layer 
assumption. The approach uses second-order spatial differencing and a 
second-order artificial dissipation model which does not require any assumption 
regarding the velocity field associated with the vortex. The equations are: 

Continuity 


dp_ 

at 


+ 


v»p u 


0 


( 1 ) 


Momentum 


d pJJ zx zi 

+ V* (pUU) = -VP + V. (x + x T ) 

8 1 


Energy 


8 ph 

JtT 


+ (puh) = -V# (Q + Q^) + — + $ + pi 


( 2 ) 


( 3 ) 


=* 

where p is density, U is velocity, P is pressure, x is the molecular stress 

tensor, x^ is the turbulent stress tensor, h is enthalpy, Q is the mean heat flux 
vector, Q'*- is the turbulent heat flux vector, $ is the mean flow dissipation rate 
and e is the turbulence energy dissipation rate. 
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In regard to the numerical method, the basic scheme used is a Linearized Block 
Implicit (LBI) ADI procedure of Briley and McDonald (1977, 1980). The splitting 
error and the linearization error associated with this basic scheme are removed by 
introducing an inner ADI iterative procedure at each time step. In brief, the 
procedure adds terms dp /dr , dpU/dr and dph/dr to the continuity, momentum and 
energy equations, Eqs. (1) - (3), respectively. Three point central differences 
are used to represent spatial derivatives and, unless specified otherwise, three 
point backward time differencing is used to represent temporal differences leading 
to a second order temporal technique. The spatial accuracy is second order except 
for the use of numerical dissipation which is discussed later. The equations, 
including the d( )/dr terms are solved by the LBI technique in a manner analogous 
to that used by Rai (1987). On convergence of the inner iteration the scheme 
becomes a fully implicit nonlinear backward time difference scheme. Extensive 
tests have been made on the most effective choice of inner iterations and physical 
time steps and these tests are discussed subsequently. 

Artificial Dissipation and Turbulence Model 

When calculating high Reynolds number flows using centered spatial 
differencing, some artificial dissipation is needed to maintain numerical 
stability and to suppress spurious oscillations in the numerical results. 

Such "artificial dissipation" could be added via the spatial differencing 
formulation (e.g., one-sided difference approximations for first derivatives) or 
by explicitly adding an additional dissipative term. The present effort utilizes 
the latter approach. When an additional term is explicitly added, the physical 
approximation being made is usually clearer than when dissipative mechanisms are 
contained within numerical truncation errors, and further, explicit addition of an 
artificial dissipation term allows greater control over the amount of nonphysical 
dissipation being added. Obviously, the most desirable technique would add only 
enough dissipation to suppress oscillations without deteriorating solution 
accuracy. Various methods of adding artificial dissipation were investigated by 
Shamroth et al. (1982), and these were evaluated in the context of a model 
one -dimensional problem containing a shock with a known analytic solution 
(one -dimensional flow with heat transfer) . The methods considered included 
second-order dissipation, fourth-order dissipation and pressure dissipation 
techniques . 
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As a result of this investigation, it^ was concluded that for nearly normal 
shocks such as that expected in this study a second-order anisotropic artificial 
dissipation formulation suppressed spatial oscillations while maintaining solution 
accuracy, and could be used to capture shocks successfully. In this formulation, 
the terms 


d — (p*- 1 d 


ax 


d<f> 

ax. 


are added to the governing equations where 4> = \x t for the x i -momentum equation and 
4> = p for the continuity equation, respectively. The exponent, n, is zero for the 
continuity equation and unity for the momentum equations. The dissipation 
coefficient, dx i? is determined as follows: the general equation without 

addition of any artificial dissipation has an x^direction convective term of the 
form ad<(>/dx i and an x.-direction diffusion term of the form 3 (ba<^/ax i )/3x^ . The 
diffusive term is expanded 


a fb3*l, a 

a — - - /axj = - — 
laxjJ 1 ax-L 


<t> ab a<p 

2 axi axi 


(4) 


and then a local cell Reynolds number, , is defined for the x^-direction by 

ab 


Re 


Ax-i 


dxi 


AXj/b 


(5) 


where b is the total effective viscosity, including both laminar and turbulent 
contributions, and Ax^ is the grid spacing. The numerical dissipation coeffic- 
ient, , is nonnegative and is chosen as the larger of zero and the local 

quantity 

^e(°x i Re Ax“ 1 ) ( 6 ) 


The dissipation parameter o xi , is a specified constant and represents the inverse 
of the cell Reynolds number below which no artificial dissipation is added. The 
maximum value of the dissipation, o Xi = 0.5, corresponds approximately to that 
obtained by using one-sided differencing throughout. A term of the form 
a(p n - 1 dx i a^/ax i )/ax i is added to each equation. 

The present results have been obtained with a mixing length turbulence model, 
and an isotropic eddy viscosity. This model assumes the existence of a mixing 


6 



length, i, and then relates an eddy viscosity, to the mixing length by 




axj ax^ J 


auj ]h 

axj J 


( 7 ) 


For flow regions upstream of the trailing edge where the flow is attached, the 
mixing length is determined by the usual boundary layer formulation 


a = /cyD i < 

•* max 


( 8 ) 


where k is the von-Karman constant, D is a sublayer damping factor and J& max is 
taken as 0.09 S , where S is the boundary layer thickness. The van Driest damping 
factor is defined as 


D = (1 - e~y + /26) 


( 9 ) 


where y"*" — y\l T /v is the dimensionless coordinate normal to the wall. 

When the mixing length formulation is used in a boundary layer environment, S 
is usually taken as the location where u/u e = 0.99. However, this definition 
assumes the existence of an outer portion of the flow where u e is independent of 
distance from the wall, and assumes that the location where u e becomes independent 
of distance from the wall marks the end of the viscous region. In an airfoil 
Navier-Stokes calculation, no such clear flow division occurs as even in the 
nominally inviscid region there are velocity gradients due to potential flow 
effects. Therefore, the boundary layer thickness, 5, is set by first determining 
u max* th e maximum velocity at each given streamwise station, and then setting S 
using 


6 


2 . 0y 


( U/U ma* = k l ) 


( 10 ) 


i.e., S is taken empirically as twice the distance from the wall to the location 
where u/u max = . In the present calculations, k t was set to 0.90. 

This boundary layer thickness specification is not suitable in the airfoil 
leading edge stagnation regions or near the rounded trailing edges where two small 
separation regions form, hence special treatment is required in these areas. The 
points 1 and 2 shown in region I of Fig. 1 are at a specified arc-length (ASj and 
As 2 ) from the nominal stagnation point, and in the present calculation ASj and 
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As 2 were approximately 5% of axial chord. In region I, a value of 6 is defined 
using a quadratic fit between the computed values and S 2 and a specified 

minimum value S min assumed at the stagnation point. The points 3 and 4 shown in 

* 

region II of Fig. 1 are at a specified arc-length (As 3 and As 4 ) from the trailing 
edge point, and in the present calculations As 3 and As 4 were approximately 6% of 
axial chord. Downstream of points 3 and 4, the value of 6 is set to S 3 and S 4 , 
respectively. In the wake region an assumed spreading rate is used. 

Obviously, both artificial dissipation and turbulence viscosity are sources of 
false diffusion and distortion of a vortex convected in a finite -difference grid. 

In the present work, the effects of artificial dissipation on the vortex structure 
is minimized by specifying a proper value of the adjustable parameter (o x ^) which 
controls the amount of added dissipation, and for the present BVI study this value 
was determined from a separate set of calculations in which the effects of this 
parameter's magnitude on the preservation of a free vortex convected over a long 
distance are examined. To minimize the spurious diffusion of vortex structure due 
to the turbulence model, the turbulent viscosity is set to be zero in a region 
ranging between the inflow section to a section which is approximately one chord 
length upstream of the leading edge of the blade. The turbulent viscosity then 
gradually blends into the values provided by the employed mixing length model. 

The numerical results obtained show these simple techniques to be effective in 
preserving the vortex during its convection towards the blade. 

Grid Construction 

During the course of this effort calculations were run with several types of 
grids. The BVI calculations require a grid which extends far upstream (perhaps 
five chords) to a station where the vortex is placed initially. Furthermore, it 
is important that the vortex be convected to the blade with minimum distortion or 
dissipation due to numerical effects. This demand argues strongly for an 
orthogonal grid upstream of the blade since significant non- orthogonality leads to 
a higher degree of numerical error and is expected to distort and dissipate the 
vortex in a physically unrealistic manner. With these considerations in mind, SRA 
has developed a "four corner" grid, a sketch of which is shown in Fig. 2. This 
grid was generated using the EAGLE code of Thompson (1987). This grid can be 
easily extended to far upstream without introducing excessive grid 
non- orthogonality in the upstream region. On the other hand, if a "C"-grid is 
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used such an extension to the far upstream would present considerable grid 
non- orthogonality in the far upstream region. Therefore a "four corner" grid is 
the grid of choice for the desired simulation, and is used here. 

As can be seen in Fig. 2, the points A, B, C, and D are points with singular 
metric information. However, they are on the boundary and the change in 

o 

coordinate line direction at these points can be kept to 45 as compared to a 
change of 90° for a standard "H"-grid. This decrease in coordinate direction rate 
of change should make the four corner grid less susceptible to numerical error in 
the vicinity of the singular boundary points. 

Figure 3 presents a computer generated plot of the overall grid for a NACA0012 

airfoil. Detailed grid distributions in the leading edge and trailing edge 

regions are shown in Figs. 4 and 5. A total of 12,000 grid points are used with 

high resolution both near the blade and in the vicinity of the blade leading 

and trailing edges. Typical resolution at the blade surface is 3.0 x 10*^ chords. 

As previously discussed, the major disadvantage with the conventional "H"-type 

grid is the presence of the singular points in the leading edge and the trailing 

edge, where the branch cut lines meet the airfoil body and the grid lines change 

their slopes at very sharp (almost ninety degrees) angles. The severity of this 

problem is reduced for the "four corner" grid by changing the location of the 

points at which branch cut lines meet the airfoil body to positions away from the 

° 

center line so that these lines can connect to airfoil with a 45 angle, as shown 
in Fig. 4. 

A calculation was performed for the isolated airfoil given in Fig. 3 at zero 
incidence. The free stream Mach number and the Reynolds number used for this 
calculation were 0.3 and 1.0x10^, respectively. The calculation was continued 
until an asymptotic steady solution was obtained. The surface pressure 
distribution at such a steady state is shown in Fig. 6 which qualitatively 
demonstrates the typical behavior of the NACA 0012 airfoil. The pressure 
distribution, for the same type of airfoil, reported by Mehta (1977) also is given 
in Fig. 6 for comparison. The results are in reasonable agreement with some 
moderate difference in the suction peak which may be attributable to the four 
corner grid. The contours of streamwise velocity, u, and transverse velocity, v, 
are given in Figs. 7-10. The major item to be noted in Figs. 7-10 are the smooth 
dependent variable fields which occur even in the presence of the coordinate 
singularities. This computation is a demonstration case for an isolated airfoil 
with a four corner grid. Comparisons against data for a variety of cases with an 
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earlier version of the code using a C-grid are given by Shamroth (1979 and 1985). 
In general, the comparisons both for pressure distribution and velocity 
distribution in the boundary layer showed good agreement with data. 

Boundary Conditions 

The authors' experience in solving Navier- Stokes equations has indicated the 
important role which correctly modeled boundary conditions play in determining 
accurate solutions and rapid numerical convergence. The present approach follows 
that of Briley, Buggeln and McDonald (1985). In brief, this approach sets total 
pressure and flow angle at the upstream inflow boundary, boundary EH of Fig. 2 and 
static pressure at the downstream outflow boundary, boundary HILE of Fig. 2. 

These represent the physical boundary conditions for the governing set of 
differential equations. In addition, the second derivative of pressure is set to 
zero at the upstream boundary, and second derivatives of all velocity components 
are set to zero at the downstream boundary. On the airfoil surface, no -slip 
conditions are applied. These are applied with sublayer resolution; wall 
functions are not used, as the first grid point off the wall is within the 
sublayer. In addition, the wall normal pressure gradient is set to zero. These 
conditions suffice for two-dimensional flow. For three-dimensional flow, as will 
be discussed later, symmetry conditions were set on the spanwise boundaries. 


III. ISOLATED VORTEX STUDY 

The first study of the present effort is an isolated vortex study aimed at 
developing a computational protocol which would preserve the vortex strength and 
shape as it is convected downstream. The present study uses three -point centered 
spatial differences as opposed to five -point upwind differences which is 
potentially more efficient and uses a numerical dissipation scheme which does not 
require any assumptions on the location or distribution of the velocity field 
associated with the vortex. However, prior to performing either two-dimensional 
or three-dimensional blade vortex interaction studies it is necessary to determine 
a run protocol which will preserve the vortex. If a numerical procedure does not 
preserve an isolated vortex, this same procedure can not be used for the blade 
vortex interaction study since non-physical dissipation will occur prior to the 
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interaction processand the process will be based upon a vortex of artificially low 
strength and large extent. 

The vortex preservation test for an isolated vortex was used by Rai (1987) as 
a measure to assess the time accuracy of different numerical schemes. In this 
test, the flow associated with a Lamb -type vortex convecting in a frees tream is 
calculated. The pressure field associated with such a vortex is a minimum at the 
vortex center and increases asymptotically to the freestream value with distance 
from the vortex center. A dissipative scheme is incapable of maintaining the 
minimum value of the pressure at the core of the vortex at its original value; 
instead, the pressure at the center increases continually as the vortex convects 
with the flow. In the absence of any physical viscosity (Euler equations) this 
change in pressure at the center of the vortex represents numerically induced 
vortex decay. Therefore, a good measure of the vortex preserving capability is 
the core pressure. Rai reported the results of a series of calculations using a 
conventional second- order- accurate , central -difference scheme (the Beam- Warming 
scheme) and an upwind biased, fifth-order scheme. Variation of vortex core 
pressure as a function of the number of core radii traveled by a vortex, as 
reported by Rai, is given in Fig. 11. Although the grid spacing in the vortex 
path was given, the overall number of grid points and the value of CFL used for 
these calculations are not reported. Under the present program a similar approach 
of monitoring the vortex core pressure as the vortex convects downstream was 
adopted to measure the time accuracy of the SRA Navier-Stokes code. A series of 
calculations were performed to assess the vortex preservation capabilities of the 
current numerical procedure and to indicate areas where refinement would be 
required. Based on the results of these calculations, the code was further 
developed to enhance the features required for the blade vortex interaction (BVI) 
study. The BVI investigation may require considerable computer resources. This 
is due to the large physical domain needed to avoid boundary effects and to the 
very fine mesh required to resolve flow details. Calculations were performed to 
study the effects of the physical boundary distances, grid spacing, artificial 
viscosity, initial vortex strength, boundary conditions and time iterations to 
optimize the use of the computer resources as a prelude to the actual BVI 
simulations . 
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Computational Grids 


Three major different types of grids were used for the free vortex 
calculations : 


1) Type "A" grid - is a rectangular grid that is 52.5 x 7.5 radii in a length to 
width with equal spacing in the x- and y-directions . The physical domain and 
the trajectory of the vortex path from its initiation to 45 radii of travel 
are shown in Fig. 12a. The grid spacing for this grid is Ax = Ay = 1/4, the 
core radius of vortex being 1.0. A portion of the grid is shown in Fig. 12b. 

2) Type "B" grid - is identical to the type "A” grid except for spacing. The 
spacing used for this grid is Ax = Ay = 1/8. The "A" and M B M grids will be 
used to demonstrate the effect of the grid density on the vortex preservation 
characteristics of the scheme. 


3) Type "C M Grid - is shown in Fig. 13. Equal spacing (Ax = y = 1/4) is used in 
the central region containing the vortex path and a stretched grid is used to 
extend the physical domain to the far field. The boundaries of the grid are 
245 x 200 radii in length to width, the core radius of the vortex being 1.0. 
The boundaries of the equally spaced grid region, the overall boundaries and 
the vortex path are shown in Fig. 14. In contrast to grid types "A" and M B" , 
the boundaries of the type "C" grid are considerably farther from the vortex 
path than the type M A" or "B" grid; calculations performed on this grid are 
intended to demonstrate the effect of the boundary extents. 


Initial Conditions. Calculation Procedures 


Initializations used for all the calculations are the same. With respect to a 
coordinate system fixed at the vortex center, the cylindrical velocity component 
of the vortex used is given by: 


v e = 


r 

2?rr 



( 11 ) 


The pressure and density for the vortex flow are obtained as below: 
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The r momentum equation gives 


. 2 
dp = pv e 

dr r 


The energy equation yields 



From Eq. (12) and Eq. (13) 



dr p 


1 

' 2 ' 

P® 7-1 Vq 

.Poo 7 2 _ 


0 


( 12 ) 


(13) 


(14) 


Eq. (14) is a first order ordinary differential equation. From Eqs . (11) and 
(14), the analytical solution of p is 

p = Pa>exp(f ) (15) 


where f is 


f - 20 

tan" 1 

f2r 2 + B 

n ' 



o 

J 4a 2 - B 2 ' 


Ll 4a 2 -B 2 ' J 



where 


D = 


1 

2 



p U 2 

00 oo 

Poo 


and 


(16) 


(17) 


B = 2a - ^ D (18) 

7 

The temperature is then given by the equation of state. It is noted here that the 
exact solution of the Lamb -type vortex satisfies the Euler equations. Options 
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were added to the code to solve the Euler equations rather than the Navier-Stokes 
equations . 

Having an analytical solution for vq , p, p and T, the initial flow field is 
established by adding a uniform flow to the above exact solution. Pressure, 
vorticity and velocity contours for the initial conditions are given in Figs. 
15-19. The flow is then allowed to develop in physical time. Since the vortex is 
translated with the freestream velocity, u^ , the location of the vortex is known 
at all times and the boundary condition is accordingly updated to the exact 
solution at each time step. Non-exact boundary conditions were also utilized in 
some of the calculations to assess the appropriate choice of boundary conditions 
for BVI problems. In these runs, the total temperature/total pressure/flow 
direction inflow boundary conditions and static pressure outflow boundary 
conditions were utilized. Most cases were run with three inner iterations per 
time step. However, tests were also conducted with two iterations per time step 
which, for the same conditions, gave identical results as those having three 
iterations. The number of iterations at each time step and the value for the 
physical time step were chosen based on the analysis performed to establish 
criteria for choosing these values. The criteria is described in later sections. 
The vortex was allowed to travel 45 radii downstream and the variation of the 
vortex core pressure was registered during the travel. At the end of vortex 
travel, plot files were written to produce contour plots. 

Central Difference Scheme. First-Order in Time 

Case #1 

A calculation was performed to determine the flow characteristics associated 
with a Lamb- type vortex convecting in a freestream, using the SRA 

first-order-in-time, central difference scheme without any inner iteration at each 
time step. The purpose of these free vortex calculations is to assess the vortex 
preservation capability of the code when applied to blade-vortex interaction 
problems. A type "B” grid was used for this case. As described before, the exact 
solution was used to set up the initial condition and the subsequent unsteady 
boundary conditions. Figure 20 shows the variation of the vortex core pressure 
with the travelled distance of the vortex. From this result it is evident that 
the first-order- time-accurate scheme is very dissipative; an improvement is needed 
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for vortex preservation. This is achieved via the second-order unsteady iterative 
implicit^ scheme described in the next section. However, the results of the Case 1 
calculation are qualitatively consistent with those obtained by Rai using the 
first-order temporal accuracy Beam and Warming procedure, although the present 

results show somewhat less dissipation. It is interesting to note that although 

» 

the vortex has lost a significant amount of its strength, the vortex structure is 
very well preserved after 45 radii of travel. This is evident from the contour 
plots of Figs. 21-24, which are taken at the final point of the vortex travel; 
i.e. after 45 core radii of travel. 

Second-Order in Time and Iteratively Implicit Scheme 

This procedure utilizes second-order accurate time-differencing, together with 
inner iterations at each physical time step. The iteration at each time-step 
decreases numerical errors resulting from ADI splitting, time linearization and 
implementation of the intermediate boundary conditions. A series of calculations 
were performed to determine the most appropriate choices for the parameters 
controlling the inner iteration pseudo- time step, the physical time step and the 
number of inner iterations. All the calculations were performed on an equally 
spaced grid distribution. The problem chosen for these tests was the 
two-dimensional vortex in a uniform flow. This is a very relevant case for the 
current overall effort and results obtained for this case are expected to have 
strong relevance for the blade vortex interaction problem. The effects of these 
new parameters on the solution behavior have been studied and are described in the 
following. 

(A) Pseudo Time-Step Parameter , DTAU 

The effect of this parameter upon the convergence behavior is shown in Fig. 

25. In this figure, the maximum residual of the streamwise momentum equation is 
chosen to assess the degree of iteration convergence. The residual of each 
equation is obtained by summing all the terms in the equation with the exception 
of the pseudo time -derivative term; i.e,, the d( )/dr term. Obviously, when the 
residual is zero, a converged solution including unsteady effects is attained. 

Since the rate of convergence depends upon the magnitude of the iteration- step , 
termed here as the pseudo time -step DTAU, it is necessary to study this dependency 
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so as to develop an efficient running, time accurate code for application to 
complex two-dimensional or three-dimensional flow cases. 

Figure 25 shows the dependence of convergence rate for the iteration at each 
time -step upon the pseudo time parameter DTAU. The calculations were initiated 
from an analytic solution and the results of Fig. 25 pertain to the first 
time-step. As can be seen, once a value of DTAU = 2.0 is reached, further 
increase does not significantly improve the rate of convergence. 

(B) Physical Time -Step DT 

After the determination of an appropriate iteration pseudo time -step DTAU, a 
physical time -step needs to be chosen to minimize the temporal error. The 
preservation of the vortex core pressure against different values of physical 
time-step (or in terms of CFL number) was chosen to determine the time accuracy 
associated with the size of the time-step. It is noted here that a time accurate 
solution should preserve the magnitude of the vortex core pressure as it is 
convected toward downstream. Obviously, numerical error will lead to some 
dissipative effects; however, it is expected that these effects can be made small. 
A series of calculations using different values of CFL were performed to determine 
this relationship. The calculations were run for the same total physical time 
interval. These runs were made with the three time-level procedure. Figure 26 
demonstrates this relationship and can be used to give guidance in choosing an 
appropriate value for the time -step for a given grid geometry. 

(C) Number of Pseudo Time Iterations 

For a specified pseudo time-step and a specified value for physical time-step, 
it is necessary to determine the number of iterations required at each physical 
time-step to achieve the desired degree of convergence. This relationship is 
shown in Fig. 27. Convergence is defined here as a state at which the per cent 
change of each and every dependent variable does not vary, at least for six 
significant digits, and the maximum residual of each and every equation reduces by 
at least four orders of magnitude during the iteration. It is noted that larger 
values of DT, the time step, not only require more iterations to achieve 
convergence (Fig. 27), they also introduce larger time truncation error (Fig. 26), 
as is expected. Using the iteratively implicit scheme, various calculations were 
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performed to assess the vortex preservation capabilities of the newly improved 
code and to study the effects of grid spacing, artificial dissipation, initial 
vortex strength and the boundary condition. 

Spatial Spacing - Two calculations were performed on grids with different spacing 
to study the effects of grid spacing on the vortex preservation. A calculation 
was also performed to investigate the effect of the physical time step. 

Case #2 

The type "A" grid was used for this calculation. Variation of vortex core 
pressure versus number of core radii travelled is shown in Fig. 28. It can be 
seen that the numerical solution is oscillating, particularly toward the later 
parts of the vortex travel. Contour plots of pressure after 45 radii of vortex 
travel are given in Fig. 29, which is a picture of a badly deformed vortex. As 
mentioned before, spacing used for the grid is Ax = Ay = 1/4. 

Case #3 

A type "B" grid was used for this calculation. The type "B" grid is identical 
to the type "A" grid except its spacing is Ax = Ay = 1/8. All the other 
conditions of the calculation were identical to Case #2. The variation of vortex 
core pressure with number of core radii is shown in Fig. 28, which indicates a 
very stable and accurate solution. It should be noted that there is almost no 
perceptible rise in core pressure after 45 radii of travel in this case. Contour 
plots of pressure, vorticity magnitude, streamwise velocity and transverse 
velocity for this case after 45 radii of travel are shown in Figs. 30-34. In 
contrast to the previous case, Ax = Ay = 1/4, this case shows no significant 
distortion of the vortex during this travel. 

Case #4 

Conditions for this case are identical to Case 3 but physical time-step At 
was doubled, i.e., it was increased from At = 0.02 to At = 0.04. This will 
correspond to increasing CFL from CFL =0.16 to CFL =0.32. The calculation was 
performed as an additional confirmation of the results presented in Fig. 26, which 
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demonstrates preservation of the vortex core pressure for different values of CFL. 
Figure 26 indicates that for CFL =0.32 the vortex core pressure is still well 
preserved. A comparison of the variation of the vortex core pressure with number 
of core radii travelled for Case 3 and Case 4 is given in Fig. 35. Contour plots 
of pressure, etc. for Case 4 were nearly identical to those of Case 3. 

Boundary Conditions : Calculations were performed to investigate the effect of the 

boundary conditions and the influence of the extents (i.e. the physical distance 
between the boundary and the vortex) of the boundaries on the vortex behavior. In 
most cases function boundary conditions were applied, i.e. the exact solution 
corresponding to the vortex moving in a freestream was imposed on the boundaries. 
This is possible since the vortex is translated with the freestream velocity, u^, , 
therefore location of the vortex is known at all times and the boundary condition 
is accordingly applied and updated to the exact solution at each time step. 

The function boundary condition discussed in the preceding section could only 
be applied to free vortex calculations. These boundary conditions cannot be 
applied to the flow calculations involving a blade. "Non- function" boundary 
conditions were used in some of the free vortex calculations as an attempt to find 
appropriate boundary conditions for the flow situations containing vortex and a 
blade. The "non- function" boundary conditions used for free vortex calculations 
are such that the total pressure, total enthalpy and flow angle are specified, and 
second derivative of pressure set to zero. It should be noted that these are time 
dependent unsteady boundary conditions and their values are changed and updated at 
each time step. Total pressure is obtained from 

p 0 = P(T 0 /T) 7/(7_1) (19) 

where 

Pq = total pressure 
P = static pressure 
Tq = total temperature 

and H 0 , the total enthalpy, is obtained from 

H 0 = h + ! 0 2 + v2 ) (20) 
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where 


h = C p T (21) 

The specifications of the above boundary conditions require the values of 
pressure p, density p , streamwise velocity u and transverse velocity v. These 
values are supplied from the exact solution. It should be pointed out here that 
the specifications of total enthalpy, total pressure and flow angle amount to 
requiring the dependent variables (i.e., u,^ v, p and T) to satisfy a system of 
algebraic equations; this is different from directly imposing functional values of 
these variables, as in the case of function boundary condition. 

The boundary conditions used on the downstream boundary are such that the 
time -dependent static pressure is specified, the second derivative of streamwise 
velocity u, transverse velocity v and the total enthalpy Hq set to zero. Again 
downstream static pressure is supplied from the exact solution. The 
"non- function” boundary conditions used on the upper and lower boundaries of the 
rectangular grid are identical to those applied on the downstream boundary. 

Case #5 

This calculation is identical to Case 2 except the type "C" grid is used 
instead of the type "A" grid. The boundaries of the type "C" grid are much 
further away from the vortex path than those of the type "A" grid (see Figs. 12a 
and 14). The results were identical to those of Case 2, indicating that when 
function boundary conditions are used, the physical locations of the boundaries 
from the vortex have little bearing on the calculated solutions. 

Case #6 

All the conditions of this calculation were identical to Case 3 with the 
exception that the alternative boundary conditions were used. Variations of 
vortex core pressure versus number of core radii is almost identical to the Case 3 
calculation (Fig. 36). Contours of pressure, vorticity magnitude, streamwise 
velocity and transverse velocity for Case 6 are shown in Figs. 37-40, 
respectively. Pressure contours of this case (Fig. 37) are identical to those of 
Case 3. However, an examination of the vorticity magnitude contours (Fig. 38) and 
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velocity component contours (Figs. 39-40) reveals the influence of the boundary 
condition on the unsteady vortex flow behavior. Compare these figures with those 
of Case 3 calculations (i.e., Figs. 31-34) where the exact function boundary 
conditions were used. It is clear that use of the alternate boundary conditions 
do distort the vorticity and velocity field in those regions where the field is 
relatively weak. For example, as shown in Fig. 39, the vorticity distortion 
occurs primarily in the region where wA> max < 0.15. The region where the major 
vortex strength occurs stays relatively undistorted. The influence of the 
boundary conditions on the vortex flow behavior can be definitely minimized by 
extending the boundaries further away from the vortex center, as they would be 
with the stretched grids contemplated for the BVI study. It should be noted that 
the upper and lower boundaries of the domain for Case 6 were only 3.75 radii away 
from the vortex center during the entire calculation, where the vortex core radius 
is 1.0. 

Initial Vortex Strength 

Two cases were conducted to assess the effect of the initial vortex strength 
upon the decay rate . 

Cases 7-8: 

The type "C" grid was used for these calculations. The initial core pressure 
of the vortex was 0.95. In all the previous cases the initial core pressure was 
0.84. In Case 7 calculations, a central difference scheme with first-order 
accuracy in time was used, whereas a second-order accurate in time and iteratively 
implicit scheme was used for Case 8 calculations. The variation of the vortex 
core pressure with number of radii travelled by vortex for Cases 7 and 8 is shown 
in Fig. 41 and for stronger vortex, having a core pressure of 0.84, in Fig. 42. 

It is clear from these figures that the weaker vortex loses less of its initial 
strength, i.e., the weaker vortex numerically preserves itself better. This is 
consistent with the fact that the degree of preservation depends upon the 
truncation error associated with difference schemes, and these errors are 
proportional to the flow gradients, which have larger values for stronger vortex. 
Pressure contours for Cases 7 and 8 after 45 radii of vortex travel are shown in 
Figs. 43-44, respectively. 


20 


Artificial Dissipation - Since all the calculations are at high Reynolds numbers, 
it is generally necessary to add "artificial dissipation" terms to suppress 
central difference spatial oscillations. Such "artificial dissipation" could be 
added via the spatial differencing formulation (e.g., one-sided difference 
approximations for first derivatives) or by explicitly adding an additional 
dissipative term. In the numerical scheme used in this study, the latter approach 
was adopted. When an additional term is explicitly added, the physical 
approximation being made is usually clearer than when dissipative mechanisms are 
contained within numerical truncation errors. Further, explicit addition of an 
artificial dissipation term allows greater control over the amount of non-physical 
dissipation being added. Obviously, the most desirable technique would add only 
enough dissipation to suppress oscillations without deteriorating solution 
accuracy. Four cases were run to assess the effect of the "artificial 
dissipatirn" on the unsteady vortex flow. In all cases the item AVISO is equal to 
o x ^ (See Eq. 6) which is approximately an inverse for all Reynolds numbers. 

Cases 9, 10, 11: 

All the conditions for these cases were identical to Case 3 with the exception 
of the amount of "artificial dissipation" added. The parameter "AVISC" is a 
measure of the amount of artificial dissipation added. The larger the AVISC 
value, the larger the added artificial dissipation. The values used in Cases 9-11 
are 0.0, 0.005 and 0.05. The value used in the Case 3 calculation was 0.001. 
Variation of the vortex core pressure versus number of core radii travelled by the 
vortex for different values of the "artificial dissipation" calculations is given 
in Fig. 45. It is clear from these curves that the higher the value of the 
"artificial dissipation" the more dissipative the calculation. As discussed 
before, there should be only enough dissipation to suppress oscillations without 
deteriorating solution accuracy. It should be noted that in the present 
calculation the solution remains stable without any added artificial dissipation 
throughout this calculation. Contours of pressure and vorticity magnitudes after 
45 radii of travel for calculations each having a different value of "artificial 
dissipation" are given in Figs. 31-32 (AVISC = 0.001) and Figs 46-51. These 
contour plots indicate that the vortex shape is very well preserved after 45 radii 
of travel, despite the fact that for higher values of "artificial dissipation" the 
vortex has lost more of its strength. 
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Solving Navier-Stokes Equations (Case 12): 


In all the cases discussed before, since a Lamb- type vortex was used, the 
Euler equations were solved. A calculation with an identical condition to Case 3 
was performed, except instead of Euler equations, Navier-Stokes equations were 
solved. The results were identical to Case 3. This is expected since, in free 
vortex calculations, in the absence of walls the contributions of the viscous 
terms are very small. 

Summary 

The original form of the SRA Navier-Stokes code which was used for this study 
was a central difference scheme, first-order accurate in time. It was shown to be 
very dissipative. The code was further developed to be second-order accurate in 
time and to do multiple iterations at each time step. These improvements proved 
to increase the time accuracy of the code significantly. Free vortex calculations 
demonstrating this aspect of the improvements were presented. Fig. 52 
demonstrates vortex preservation capabilities of the original SRA central 
difference scheme code and the improved second- order accurate in time and 
iteratively implicit scheme code. Results of Rai's fifth-order-accurate upwind- 
biased scheme are also presented in this figure for comparison. Numerous test 
cases were performed to assess the effects of grid size, initial vortex strength, 
boundary conditions and artificial dissipation on the unsteady free vortex flow 
calculations. The calculation conditions and the results are summarized in Table 
1. It can be noted from Fig. 52 that the inclusion of the second order accurate 
scheme, together with the use of iteratively implicit techniques greatly increase 
the time accuracy of the solution and the numerical simulation and preserve the 
vortex strength and structure for the duration needed for the blade vortex 
interaction study. 
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IV. TWO-DIMENSIONAL BLADE VORTEX INTERACTION 


The second case considered is that of a two-dimensional blade vortex 
interaction. The isolated vortex study showed the proper choice of time step, 
iteration step, iteration number, etc. required to preserve the vortex strength 
and radius. However, this study does not address resolution and artificial 
dissipation issues for the BVI problem. These issues were addressed in a 
two-dimensional BVI simulation. 

Grid. Boundary Conditions and Artificial Dissipation 

Fig. 53 illustrates the grid distribution used for the BVI simulation. The 
total number of grid points is 144 x 118. The inflow boundary is located at 7 
chords from the blade leading edge while the outflow boundary is located at 5 
chords from the blade trailing edge. Based upon the isolated vortex study the 
distance between the top boundary and the chord line of the NACA0012 airfoil was 
set at 5 chord lengths. The geometric configuration is symmetric about the chord 
line. Along the inflow boundary, the total pressure, the total temperature and 
the inflow angle are specified. The pressure is obtained by extrapolation. Along 
the outflow, top and bottom boundaries, the static pressure is specified, the 
velocity and the total temperature are obtained by extrapolation. On the blade 
surface, non-slip conditions are imposed. The density is obtained by solving the 
continuity equation and the surface temperature is specified as the constant, free 
stream total temperature. When the vortex is introduced into the transonic flow 
field at a point upstream of the blade, the boundary values must take into account 
the existence of this vortex, e.g., the inflow angle will not be zero and 
generally is not uniform along the inflow boundary. In addition, as the vortex is 
convecting towards the blade, these boundary values are changing with time. For 
the present simulation, they are taken from the composite "vortex in a 
free -stream" solution. 

A spatially varying artificial dissipation was used in the BVI calculation 
presented here. Upstream of the blade AVISO was set to 0.005. Based upon the 
vortex preservation studies this should keep the vortex decay within allowable 
levels prior to the interaction. However, it should be recognized that the vortex 
preservation studies were done on a uniform grid whereas the present calculation 
is done on a non-uniform grid. Downstream of the airfoil leading edge station the 
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value of AVISO was set at 0.05 except in the immediate vicinity of the leading and 
trailing edge locations where it was set to 0.5. The relatively high value of 0.5 
was required due to the highly stretched, highly skewed grids and the marginal 
resolution of the flow in these regions . The regions of different AVISO values 
were blended at boundaries to avoid discontinuous changes in the artificial 
dissipation. 

Although this procedure for specifying numerical dissipation does put a higher 
than desired level in the immediate vicinity of the leading and trailing edge 
regions, it is necessary to suppress oscillations in these regions. This same 
procedure has been used in a variety of airfoil and cascade studies and has given 
good agreement with experimental data for a wide variety of cases, e.g. Shamroth 
(1985), Weinberg et al . (1986) and Shamroth et al . (1988). 

Flow Parameters 

The reference length is the chord length of the blade and the reference flow 
conditions are the free stream condition with M*, — 0.8 and Re = 1.0 x 10**. The 
background flow is a steady transonic flow with shock waves standing in the middle 
of the blade. Furthermore, the flow is symmetric about the chord line; hence, the 
lift coefficient (C^) is zero. The surface pressure distribution of this 
background flow is shown in Fig. 54. 

The dimensionless strength and core radius of the vortex are —1.6 and 0.2, 
respectively, where the minus sign indicates that the vortex has a clockwise 
sense. The initial location of the vortex center is at a point 5 chords upstream 
of the airfoil leading edge (x v = —5.0) and 0.26 chords below (y v = —0.26). The 
calculation is carried out from t = 0 to t = 8 with constant time step At = 0.005. 
It is noted here that the vortex core arrives at the blade leading edge when 
t = 4.95, which indicates an average core velocity of 0.99 v^. 

Results and Discussion 


The time histories of the aerodynamic coefficients during the blade-vortex- 
shock interactions are described in Figs. 55-57. It should be noted that the 
small amplitude variations occurring for t up to 1.0 are due to the impulsive 
introduction of the vortex into the background flow, and these small amplitude 
oscillations have been damped out long before the onset of significant blade 
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vortex interactions. Furthermore, these coefficients are evaluated in terms of 
static pressure; they do not include the contribution of viscous stresses. The 
lift coefficient (C L ) and the quarter-chord pitching moment coefficient (C^) are 
shown in Figs. 55 and 56, respectively. It should be noted that the background 
flow is a nonlifting case and that any lift generated during the interaction is 
induced by the vortex. Since the vortex flow is revolving in the clockwise 
direction, when the vortex is approaching the blade, it induces nonuniform and 
unsteady velocities that result in negative angles of attack at the blade. This 
influence changes to increasing angles of attack after the vortex has reached the 
blade. Severe load variations occur during the time period from t = 4.0 to 
t = 6.0; i.e., when the vortex is within one chord length of the blade leading 
edge. During this period of time, and change their signs, while C D 
undergoes rapid variation exhibiting two distinct temporal maxima (Fig. 57). The 
coefficients have not yet returned to their undisturbed value at t = 8.0 
indicating that the interaction is still affecting the flow. This is consistent 
with the results of Srinivasan et al. (1986). 

The interaction between the vortex and the blade with a shock are further 
elucidated in terms of the instantaneous static pressure distribution at several 
selected time stations. Fig. 58(a) gives the pressure contours over the entire 
computational domain at t = 0, while Fig. 58(b) gives the distribution of static 
pressure coefficient on the blade surface at t = 0. This is the starting flow 
field. As the vortex convects towards the blade, the upper surface shock moves 
in the upstream direction, its strength is decreasing and the extent of the 
associated supersonic pocket also is reducing. On the other hand, the lower 
surface shock moves in the downstream direction with increased strength. In 
addition to the motion of the shock waves, pressure difference between the upper 
and lower surfaces start to build up. These generic features are illustrated in 
Fig. 59(a) and (b) , which are obtained at t = 2.0. The outer pressure contours on 
the aft region of the blade shows some "wiggles"; however, this is in a region of 
nearly uniform pressure and, therefore, the presence of these "wiggles" is not 
significant, but represents small changes in a nearly uniform field. The flow 
field at t = 4.0, i.e., when the vortex core is about one chord length upstream of 
the blade, is shown in Fig. 60(a) and (b) . The upper surface supersonic pocket 
practically has disappeared. The lower surface shock wave becomes stronger and is 
located in a further downstream position; at the shock's root the flow shows signs 
of separation. In addition, a significant transverse pressure gradient exists in 
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the leading edge region. At t = 4.5, this leading edge transverse pressure 
gradient becomes the dominant feature of the interaction. In addition, the lower 
surface shock stops moving towards the trailing edge, in spite of the fact that 
its strength is still increasing. These features are depicted in Figs. 61(a) and 
(b) . When t = 5.0, the vortex core 'hits' the blade, the averaged Mach number 
inside the lower surface supersonic pocket increases to approximately 2. Most of 
the disturbances on this surface will propagate downstream until reaching the 
shock. Part of the disturbance is able to leave the supersonic pocket near the 
outer region of the shock. However, the remainder seems to build up at the root 
of the shock, as indicated by Fig. 62(a). At the same time, high pressure 
disturbance starts to be released from the upper surface of the leading edge, as 
indicated by Fig. 62(a). The state of the flow at this moment is very volatile, 
within a short period of time, this process of disturbance build-up has collapsed, 
as illust ated by Fig. 63 at t = 5.5. The shock is no longer an approximately 
normal shock but has an oblique leg which intersects the blade wall upstream of 
the intersection location at t = 5.0 (Fig. 62(a)). The emission of a high 
pressure pulse from the upper surface of the leading edge is evident from Fig. 
63(a); this high pressure pulse then propagates upstream in a domain including the 
frontal region of the entire leading edge (see Fig. 64(a)). Between this frontal 
high pressure region and the lower surface shock wave, a low pressure pulse is 
propagating towards the lower outer boundary. Figs. 62(b), 63(b) and 64(b) 
illustrate the collapse of the disturbance building up process and the subsequent 
relaxation of the strength and the location of the shock wave on the lower 
surface. The general features of the flow at t = 6.0 are: the existence of a 
supersonic pocket on the lower surface, significant flow separation originating 
at the root of the shock, the appearance of vortex remnants near the blade 
trailing edge, and the development of supersonic flow on the upper surface. 
Subsequently, the lower surface supersonic pocket continually reduces its extent 
and eventually disappears by t = 8.0. The lower surface shock moves towards the 
leading edge with continually diminished strength and then vanishes. At t = 8.0, 
the vortex remnants have been further convected downstream. The flow on the lower 
surface does not exhibit any appreciable separation and is entirely transonic. 
Furthermore, about 70% of the upper surface is covered by a supersonic pocket, 
with compression waves appearing near the trailing edge of the blade. It is clear 
that the interaction is a strong one and the vortex path particularly, when the 
vortex is near the blade is influenced by the blade and the interaction. The 
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local velocity changes significantly during the interaction. This makes 
techniques based upon assumed vortex position or shape unlikely to provide an 
accurate simulation. 

The radiation of pressure pulses from the leading edge region as a result of 
the blade -vortex interaction is further investigated in terms of the scaled 
pressure disturbance. In the two-dimensional linear far field, the amplitude of 
propagating waves should be proportional to 1/r 1 *, where r is the distance to the 
origin of the disturbance. The scaled pressure disturbance is obtained by 
subtracting the steady- state pressure background value from the instantaneous BVI 
solution and then dividing the difference by the far upstream value of pressure. 
The result is then multiplied by r h , where r is the distance to the leading edge. 
Fig. 66 defines the points at which the time histories of the scaled disturbances 
are recorded and presented in Figs. 67 (a), (b) and (c) . These points are taken 
as grid points leading to some variation in d. The first peak, I, is due to the 
interaction between the incoming vortex and the leading edge, when the vortex 
center is within one chord length of the blade. The valley, II, is associated 
with the passage of the vortex core through these points, and the peak, III, 
originates from the upper surface of the leading edge, after the vortex case 
'hits' the blade, as described before. The results of Fig. 67 shows the 
simulation to give 1/r ** scaling only approximately. Points #1 and #2 satisfy the 
scaling reasonably well. The results for point #3 show the same shape but the 
level of the maximum and minimum distributions are significantly more negative 
than points 1 and 2. 

The lift coefficient time history of the present effort showed the same 
qualitative features as that obtained by Srinivasan, McCroskey and Baeder (1986) 
and Sanker and Tang (1985) and the moment coefficient showed qualitative 
similarity with that of Srinivasan et al. (1986). Surface pressure distributions 
also showed qualitative agreement with the results of Srinivasan et al. (1986) and 
Rai (1987) as the shock strength on the upper surface decreased and that on the 
lower surface increased during the interaction. Finally, the details of the 
acoustic wave leaving the blade during the interaction show considerable 
similarity with those of Rai (1987). 
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V. THREE-DIMENSIONAL BLADE VORTEX INTERACTION 

The final case considered is that of a three-dimensional blade vortex 
interaction. The two-dimensional isolated vortex study clearly indicates the need 
for appropriate grid resolution to avoid unacceptable levels of numerical 
dissipation in the vortex. This need to maintain high grid resolution even in the 
three-dimensional case influenced the choice of case to be studied. Although 
numerical simulation of a blade with a tip is the eventual goal, inclusion of a 
tip would require resolution of an entire new set of length scales associated with 
the tip geometry and the local tip flow; this would increase the number of grid 
points required beyond what is practical for this effort. Therefore a simpler 
problem which contains the basic three-dimensional flow field physics was chosen 
for the present investigation. 

A ske ch of the case to be considered is shown in Fig. 68, which shows a 
periodic vortex initially consisting of straight line segments in a plane parallel 
to the airfoil midline. Since there is a natural symmetry, the computational 
domain is reduced to the region between spaced spanwise symmetry planes and, 
further, since no tip is included an equally spaced spanwise grid can be used. 

This allows use of a computational grid consisting of 21 equally spanwise planes 
with each plane containing 144 x 118 points for a total of 356,832 grid points in 
the computational domain. The specific configuration chosen had an angle <p of 20 
degrees, a distance d of 2.0 chords and points a and b were 5.0 and 4.27 chords 
upstream of the blade, respectively. 

The equations solved are the three-dimensional Navier- Stokes equations as 
given in Section II. These solve the equations in an inertial frame. Extension 
to a rotating frame can be made either by having the wing move in this inertial 
frame as was done in the dynamic stall study of Shamroth (1985) or by adding 
centrifugal and Coriolis terms. This latter approach was used at SRA in a 
three-dimensional rotor-stator interaction study using the same computer code as 
used here by Gibeling et al. (1990). It is this latter approach which would be 
recommended if the rotational effects were to be included. Although this 
additional capability was not available at the initiation fo the present study, it 
is now available and adds little complexity to the numerical solution. Obviously 
boundary conditions must be written in the rotating frame but this should be 
straightforward. 
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Having specified the initial vortex geometry, it is necessary to create 
compatible velocity, density and pressure fields. Considering a vortex line 
segment, as shown in Fig. 69, the induced velocity at a point P due to this line 
segment is given by: 


AV' = 
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"relative" position vector of P defined from the end points to P. 
minimum distance of P from the vortex line, 
distance between the end points "1" and "2" of the line 
element. 

a constant and its sense is consistent with a vorticity vector 
directed from "1" to "2". 


The influence of the core is accounted for by multiplying AV by a factor 

„ 2 f 2 2 

f core ~ r m / ^ r m + r c 

where r c is the core radius and it is the location at which the tangential 
velocity has its maximum. Induced velocity due to the overall geometry consisting 
of many line segments is given by: 

V' = lIV' (23) 

The number of filaments is chosen such that the induced velocity field satisfies 
the symmetric condition for the flow in the computational domain. Experience 
indicates that, by placing approximately 30 or more vortices on each side of the 
representative vortex, an excellent approximation to the spanwise periodicity 
required by an infinite number of vortices can be achieved. Fig. 70 shows the 
overall geometry and the relative position of the vortex and the blade. It should 
be noted here that the induced flow field associated with the filament present in 
the computational domain contains not only the dominant contribution of this 
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vortex but also the contributions of all the other vortices in the same array. 

The vortex model constructed according to Eq. (22) and illustrated in Fig. 70 is 
to be used only in setting the initial flow condition of the simulation. The 
subsequent flow development is governed by the solution of the Navier-Stokes 
equations . 

The composite flow will consist of the background flow and the variation from 
that background flow due to the initiation of the vortex filament. The presence 
of the blade wall is incorporated into this initial flow specification by 
multiplying the induced velocity field v' by tanh (Acd/S) where: 

d is the distance to the wall 
K is the von Karman constant 

and 

S is the boundary layer thickness 

thus 

u' = u' tan h (tcd/6) (24) 


and u' satisfies the no-slip condition. 

The composite velocity field is: 

(25) 

(26) 
(27) 


U = Ujr, + ll' 

v = v b + v' 

W = Wfc + w' 


where u^, vjj and w^ are streamwise velocity, spanwise velocity and transverse 
velocity of the background flow; i.e., the flow without a vortex. 

The initial pressure field of the composite flow must be evaluated in such a 
way that it is consistent with the prescribed composite velocity field. Two 
approaches were attempted. Under the first approach the temperature was assumed 
to be that of the background flow. 


T = T b (28) 

The strategy was to find the density ( p ) and pressure (P) of the composite field 
such that they are consistent with the velocity field given by Eqs. (25) - (26) 
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and the temperature field given by Eq. (28). This was done by initially assuming 


P = Pb 

where p ^ is the background density, and then integrating the suitable component of 
the momentum equation in its steady form 


t 1 - 6 - ft “ °] 


to obtain the pressure field. Utilizing this pressure, temperature (T) and the 
equation of state, p can be updated. Using the updated p and integrating the same 
component of the momentum equation, a new pressure (P) is obtained. Convergence 
of this iteration process implies that a consistent P, p field is achieved. 
Although a convergent solution could be obtained the numerical truncation error of 
the integration process was excessive and an alternative numerical specification 
was sought . 

Under this alternative the flow field was constructed on a spanwise plane by 
spanwise plane basis. At each spanwise plane the vortex was assumed to be a Lamb 
vortex and the pressure, density and temperature found via Eqs . (11) - (18). 
Although this does represent an approximation, it did provide an initial starting 
flow field from which the computation was able to proceed without any significant 
discontinuity. 

The initial pressure field is shown in Fig. 71. A top view of the pressure 
contours is shown in the upper portion of the figure with the pressure contours 
due to the vortex being on the left hand side of the figure and that of the blade 
on the right hand side of the figure. The inclination of the vortex relative to 
the blade is evident. The lower portion of the figure shows a perspective plot 
with contours plotted in three spanwise planes. It should be noted for future 
reference that the vortex is initially closest to the blade at LY = 21, and 
furthest upstream from the blade at LY = 1. The initial pressure and Mach number 
fields through the center spanwise plane, LY = 11, are shown in Fig. 72 and Fig. 

73. There is no spanwise pressure variation for this initial condition. 

The remaining figures present either pressure contours or Mach number contours 
at a specific time in three spanwise planes or surface pressure distributions on 
the blade. When figures showing contours are vertical, the 
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LY = 1 plane which initially has the vortex at its furthest upstream location is 
at the top of the figure, LY = 11 is in the middle and LY = 21 is at the bottom. 
When presented horizontally, LY = 1 is on the left. 

Figure 74 shows pressure contours at t = 1.23, At = 1 is the time a particle 
at undisturbed frees tream velocity traverses a distance of one chord. At all 
spanwise locations the movement of the upper shock forward and the lower surface 
shock backwards is evident. The vortex centerline which was initially in a plane 
parallel to that of the blade midspan has migrated out of that plane, with the 
vortex core location being at progressively higher locations as the spanwise 
location proceeds from plane LY = 1 to plane LY = 21; i.e., the portion of the 
vortex for that downstream rises due to the self- induced motion as is expected 
(e.g. Hama and Nutant, 1963). The result is an interaction where the vortex is at 
angle to the blade and is no longer parallel to the blade centerline which 
changes vortex location relative to the blade as a function of span. 

Figures 75 and 76 present pressure and Mach number contours at t = 2.838. At 
this time the upper surface shock has nearly disappeared at. all spanwise 
locations.. The height of the vortex core clearly varies with spanwise location 
and a strong shock remains on the lower surface with the shock appearing to be 
significantly stronger than that shown at t = 0.0, Fig. 72. In all cases the 
contours are quite clear and free of oscillations. A composite surface 
pressure plot consisting of surface pressure plots at 7 different spanwise 
locations is shown in Fig. 77; no significant spanwise pressure variation is 
observed at this time. Pressure and Mach number fields at t = 4.139 are 
presented in Figs. 78 and 79. At this time the interaction is very strong at all 
streamwise locations. The interaction at LY = 21 has proceeded the furthest. A 
composite surface pressure plot is shown in Fig. 80; at this time significant 
surface pressure differences are evident with the highest suction peaks occurring 
at the LY = 21 location. This is the spanwise location at which the vortex is 
closest to the blade and at which the vortex is highest. The shock location 
varies mildly with span; however, the shock strength is approximately constant 
with span. The vortex center at the LY = 21 plane has nearly reached the blade 
leading edge at this time. The fields at t = 4.56 are presented in Figs. 81 and 
82. At t = 4.56 the three-dimensionality of the interaction is very clear from 
both pressure and Mach number plots. By t = 4.56, the shock on the lower surface 
at the LY = 21 plane has moved forward whereas that on the LY = 1 plane remains 
aft on the blade. The surface pressure distributions indicate shock bifurcation 
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at some spanwise locations. As in the two-dimensional case, high pressure pulses 
are being released from the upper surface of the leading edge at all spanwise 
locations at this time. The pulses are most apparent at the LY = 21 plane; 
however, they are apparent even at the LY = 1 plane. A composite pressure plot is 
shown in Fig. 83. 

At time, t =4.7, the interaction is very strong over the entire spanwise 
extent of the domain. The pressure contours are shown in Fig. 84 and the Mach 
number contours in Fig. 85. Based upon these figures the lower surface shock 
system on the LY = 21 plane has moved considerably forward from its most aft 
position and continues to sharpen in the freestream. The LY = 1 system is a 
single shock on the aft portion of the blade. Figs. 86-88 show details of the 
contours. These plots were created from a culled plot file hence the somewhat 
discontinuous appearance of the contours. Of particular interest is the vortex 
interaction induced shock bifurcation shown in Fig. 88. This represents the 
spanwise plane where the interaction has proceeded furthest. The contours at LY = 
11 in Fig. 87 show bifurcation starting, whereas those in Fig. 86 show no 
significant sign of bifurcation. Pressure and Mach number fields at t = 5.157 are 
given in Figs. 89 and 90. At this location the effects of three dimensionality 
are significant, as can be seen by comparing Figs. 81 and 89. Figure 81 presents 
results at t = 4.56, at which time the vortex at spanwise plane LY = 21 has just 
passed the airfoil leading edge. Figure 89 presents results at t = 5.15, at which 
time the vortex at LY = 1 has just passed the leading edge. At LY = 1, t = 5.15, 
(Fig. 89) the sfhock shows a clear bifurcation which is not present at LY = 21, t = 
4,56 (Fig. 81). This is likely to be an effect associated with the three 
dimensionality. Details of the lower surface pressure contours at t = 5.157 are 
presented in Figs. 91-93, where the different shock locations and strengths at the 
various spanwise locations are apparent. A composite surface pressure plot is 
given in Fig. 94 and lower surface contour plots are shown in Fig. 95. 

Plots of blade lift, drag and moment coefficient as a function of time are 
given in Figs. 96-98. Although these are qualitatively similar to their 
two-dimensional counterparts, they are affected by the variation of vortex height 
with spanwise location. Surface pressure contours on both the upper and lower 
surfaces are shown in Figs. 99-107. 

The general qualitative picture obtained for the three-dimensional interaction 
is similar to that previously obtained for the two-dimensional interaction as 
modified by the inclination of the vortex relative to the blade leading edge and 
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the migration of the vortex core line so that it does not lie on a plane parallel 
to the blade midplane. As the vortex approaches the blade, the core location at 
one spanwise boundary rises and that at the other spanwise boundary falls. At all 
spanwise locations the upper surface shock moves forward, its strength reduces and 
the associated supersonic pocket reduces. The extent of the interaction at any 
specified time is a definite function of the spanwise location. However, although 
at any spanwise location the general flow development is similar to that of the 
two-dimensional case; three-dimensional effects influence the flow details. The 
lower surface shock moves aft with increasing strength. As the vortex core moves 
closer to the airfoil the upper surface shock system disappears. All these 
interactions, as well as the ones discussed subsequently, occur first in the 
LY = 21 plane, which is the plane in which the vortex is initially closest to the 
blade. During the interaction the lower surface shock system reaches a maximum 
downstrerm position and then moves forward and sharpens. High pressure 
disturbances appear at the leading edge and propagate upstream. As the 
interaction continues, the shock shows a bifurcation with the appearance of an 
oblique leg, which produces a diffuse pressure rise at the blade surface; however, 
the shock remains very sharp in the outer flow removed from the blade. The 
pressure pulse continues to propagate in the upstream direction. 

The present BVI simulation is based upon a solution of the full Navier-Stokes 
equations. The only previously published viscous solution to this problem to 
the present authors' knowledge is that of Srinivasan and McCroskey (1989) based 
upon the thin shear layer approximation to the Navier-Stokes equations. The 
Srinivasan and McCroskey effort considered unsteady interaction of a rotor with a 
vortex with the equations solved in an inertial frame with the rotor being 
represented as a moving no-slip/no through flow boundary. A presecribed vortex 
formulation was used. The present approach uses the full three-dimensional 
Navier-Stokes equations without any thin shear layer assumption and does not split 
the velocity field into a prescribed vortex portion and a background flow. 

However, the present effort considers a simpler problem, that of a vortex 
encountering a blade section of finite spanwise extent in an inertial frame. 
Therefore, no direct comparisons of the results of the two efforts can be made. 

Based upon the three-dimensional simulation of the present effort certain 
conclusions can be made. For vortices of modest angle to the blade, the present 
case was 20° , the three-dimensional interaction is qualitatively similar to the 
two-dimensional interaction on a section by section basis. However, differences 
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do exist as details of the flow field at any spanwise locations for which the 
interaction had occurred at any earlier time. This effects the details of the 
shock shape, shock location, surface pressure distribution and field pressure 
distribution. Therefore, although qualitative conclusions can be drawn from 
two-dimensional simulations, a need for detail or quantitative conclusions require 
a three-dimensional simulation. This conclusion is in basic agreement with that 
of Srinvasan and McCroskey (1989). 

CONCLUDING REMARKS 

A Navier-Stokes analysis has been applied to the problem of two- and 
three-dimensional blade vortex interactions. The two-dimensional interaction 
studies have shown the ability of an iterative implicit procedure using three 
point central differences to provide a highly non-dissipative solution procedure. 
Results were obtained which showed details of the flow including the propagating 
pressure disturbances resulting from the interaction. The results showed many of 
the known physical features and gave good qualitative correspondence with other 
analyses. The three-dimensional study considered a blade vortex interaction with 

o 

the vortex initially at 20 incidence to the blade. A solution was obtained 
showing general features of the two-dimensional interaction with modifications due 
to the angle between the vortex and the blade. As in the two-dimensional case, 
pressure disturbances were observed leaving the blade leading edge. Details of 
the simulated flow physics are presented in the appropriate sections of this 
report. 
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SYMBOLS 


C D 

C L 

C M 

C P 

D 

d 

DT 

D TAU 
dx^ 

f 

'-core 

h 

H 0 

A 

LY 

P 

Q 

r 

T 

t 


u i 

u, v 

u b’ v b’ w b 


V 

v e 

x i 

y + 

r 

s 

Ax, Ay 

Ax^ 

e 

K 


Vortex core radius 

Drag coefficient 
Lift coefficient 

Moment coefficent 
Specific heat 
Damping coefficient 
Distance from wall 
Time step 

Pseudo time step parameter 
Artificial dissipation coefficient 
Velocity damping factor in vortex core 
Enthalpy 

Stagnation enthalpy 
Mixing length 
Spanwise plane number 
Pressure 

Heat flux vector 
Radial distance 
Temperature 
Time 

Velocity component in i^ direction 
Cartesian velocity component 

Background flow velocity components excluding contribution of 

vortices 

Velocity 

Tangential Velocity 
Cartesian Coordinate 

Dimensionless coordinate normal to wall 
Circulation 

Boundary layer thickness 
Grid spacing 

Grid spacing in i^ direction 
Turbulence dissipation 
Von-Karman constant 
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n 

7T 


P 

T 

$ 

<P 


Viscosity 
Stress tensor 
Density 

Relaxation time step coordinate 
Mean flow dissipation 
Dependent variable 


Subscripts 


b 

o 

00 


Background value 
Stagnation value 
Far upstream value 
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Fig. 1 Nomenclature for Leading and Trailing Edge Turbulence Length 
Scale Specification 
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Grid Distribution for NACA 0012 Airfoil 
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Fig. 4 Leading Edge Grid Distribution 




45 


Fig. 5 Trailing Edge Grid Distribution 
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Fig. 7 Streamwise Velocity Contours 
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Fig. 9 Transverse Velocity Contours 
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Fig. 12 Physical Domain and the Vortex Path for the "A" Type Grid. 
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Fig. 13 Grid Spacing for Type "C" Grid 
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Fig. 14 Physical Domain and the Vortex Path for Type "C" Grid. 
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Fifc. 15 Contours of Pressure at Initialization 
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Fifc. 16 Contours of Pressure at Initialization 
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Fig. 17 Contours of Vorticity Magnitude at Initialization 
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Fig. 19 Contours of Transverse Velocity at Initialization 
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Fig. 21 Contours of Pressure after 45 Core Radii of Travel 

(first-order-time-accurate central-difference scheme) 
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Fig. 22 Contours of Vorticity Magnitude after 45 Core Radii of Travel 
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Fig. 23 Contours of Streamwise Velocity after 45 Core Radii of Travel 
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Fig. 24 Contours of Transverse Velocity after 45 Core Radii of Travel 
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Versus DTAU After 4 Iteration 


89*0 



CD 

ID 


cn 

Csj 

CD 

CD 

CD 

CO 

CD 

• 

• 

• 

• 

• 

O 

O 

O 

o 

O 


(°°d/d) 3UflSS3dd 3dOO X31UOA 


80 


Fig. 26 Effects of Tim- 
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Fig. 27 No. of Iteratii 
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Fig. 29 Contours of Pressure after 45 Core Radii at Travel 
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Fig. 30 Contours of Pressure after 45 Core Radii of Travel 
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Fig. 31 Blow up of Contours of Pressure after 45 Radii of Travel 
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Fig. 32 Contours of Vorticity Magnitude after 45 Core Radii of Travel 
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Fig. 33 Contours of Streamwise Velocity after 45 Core Radii of Travel 
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Fig. 3A Contours of Tranverse Velocity after A5 Core Radii of Travel 
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Fig. 37 Contours of Pressure after 45 Radii of Travel - Alternate Boundary Conditions Used 
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Fig. 38 Contours of Vorticity Magnitude after 45 Core Radii of Travel - Alternate Boundary Conditions Used 
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Fig. 39 Contours of Streamvise Velocity after 45 Core Radii of Travel - Alternate Boundary Conditions Used 
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Fifc. 40 Contours of Transverse Velocity after 45 Core Radii of Travel - Alternate Boundary Conditions Used 
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Fig. A3 Pressure Contours after A5 Radii of Vortex Travel for Weaker Vortex. 
Central Difference Scheme - First Order Accurate in Time. 
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Fig. 46 Pressure Contours after 45 Core Radii of Travel, AVISC=0.00 
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Fig. 47 Vorticity Magnitude Contours after 45 Radii of Travel, AVISC=0.00 
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Fig. 48 Pressure Contours after 45 Radii of Travel, AVISC= 0.005 
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Fig. 51 Vorticity Magnitude Contours after 45 Radii of Travel, AVISC=0.05 
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Fig. 53 Grid Distribution for BVI Simulation. 
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Fig. 56 Time History of the Quarter-Chord Pitching Moment Coefficient, 
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Fig. 57 Time History of the Drag Coefficient, 



Fig. 58a Static Pressure Contours (t-0) 
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Fig. 58b Surface Cp (t«0) 
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Fig. 59a Static Pressure Contours (t=2.0) 
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Fig. 59b Surface Cp (t=2.0) 
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Fig. 60a Pressure Contours (t=4.0) 
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Fig. 60b Surface Cp (t-4.0) 
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Fig. 61a Pressure Contours (t =s 4.5) 
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Fig. 61b Surface C p (t-4.5) 
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Fig. 62a Pressure Contours (t-5.0) 
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Fig. 62b Surface Cp (t**5.0) 
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Fig. 63a Pressure Contours (t=5.5) 
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Fig. 63b Surface C p (t-5.5) 
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Fig. 64a Pressure Contours (t-=6.0) 



Fig. 64b Surface C p (t«6.0) 
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Fig. 65 Surface 
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Fig. 66 Points where Scaled Pressure Disturbance are Calculated. 
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Fig. 67a Scaled p' at Point #1. 
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Fig. 67b Scaled p' at Point #2. 
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Fif>. 67c Scaled p' at Point #3. 
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(-7.0, -5.0) (b) X-Z Plane. 

Fig. 68 Physical Domain for 3D BVI Simulation. 



Figure 69 Vortex Filament. 


144 



Figure 70 Overall Geometry. 
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Fig. 71 Initial Pressure Field, t*0.0. 
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Fig. 72a Initial Mach Number Field, t~0.0. 


CONTOUR LEUELS 
0. 3-1000 
0.36000 
0.38000 

0. -10000 
0. -12000 
0.^000 
0. 16000 
0.-18000 
0.50000 
0.52000 
0.5-1000 
0.56000 
0.58000 
0.60000 
0.62000 
0.6-1000 
0.66000 
0.68000 
0. 70000 
0. 72000 
0. 7-1000 
0. 76000 
0. 78000 
0.80000 
0.82000 
0. 8-1000 
0.86000 
0.88000 
0.90000 
0.92000 
0.9-1000 
0.96000 
0.98000 
. 00000 
. 02000 
.0-1000 
06000 
UHOOO 


1 . 



Fig. 72b Initial Pressure Field, t-0.0 
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Fig. 75 Pressure Field, t«2.838. 
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Fig. 76 Mach Number Field, t«2.838. 
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Fig. 77 Composite Surface Pressure, t-2.838. 
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Fig. 79 Mach Number Field, t«4.138. 
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Fig. 80 Composite Surface Pressure, t-4.138. 
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Fig. 82 


Mach Number Contours, t-4.56. 
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Fig. 83 Composite Surface Pressure, t-4.56. 
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Mach Number Field, t-4.7. 
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Fig. 87 Details of Suction Side Pressure, LY-11 
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Fig. 89 Pressure Field, t-5.157 
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Fig. 90 Mach Number Field, t-5.157. 
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rig. 92 Details of Suction Sid® Pressure, LY-ll, t-5.157. 
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Fig. 93 Details of Suction Side Pressure, LY-21, t-5.157 
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Composite Surface Pressure, t-5.15 
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Fig. 95 Pressure Field at t = 5.157 
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Fig. 96 Blade Lift Coefficient. 


o 



ooooooooo 

1N3I 31 33303 9byQ 


4 -> 

c 

<D 


e'- 

er* 

Ui 


198 


PHYSICAL TIME CT1 


0.40 


o 



oo cm o ^ o\i on <^r 

«••••••• 

oooooooo 

I I I I 

1N3I 31 33303 1N3W0W 



<x 

(_) 

cn 

>- 

x 

<L 


199 


Fig. 98 Blade Moment Coefficient 
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Fig. 99 Surface Pressure Contour at 
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Fig. 100 Lower Surface Pressure Contour, t-4.138. 


CONTOUR LEUELS 



c 

c 

B 

B 

B 

B 

C 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 


E 

E 

B 


B 

B 

B 

B 

B 

B 

B l 


B 

B 

B 

B i 


B 

B 

B 

3 

■ 

X 

B 

X 

X 

X 

X 

B 

x 

~ 

x 

B 

X 

X 

X 

X 


B 

B 

X 

E 

E 

E 

E 

E 

E 

E 

E 


E 

E 

E 



E 

E 

x 

x 

x 


- 

X 

x 

B 

X 

B 

x 

— 

B 

6 

X 

X 

X 

B 


B 

B 

E 

E 

E 

E 

E 

E 

E 

E 

E 


E 

E 

E 



E 

E 


— 

2 


B 


T 

X 

E 

B 

Ps. 


X 

ec 

C 

N 


X 

X 

X 

cs 


X 

E 

E 

PN. 

- 

X 

E 

E 


1 

X 

E 


Ps 

T 

X 

o 

B 

B 







esi 

C\ 

esi 

N 

N 

r. 

n 

n 

n 

r. 


T 


T 

*T 

U! 

in 

in 

m 

in 

X i 

X 

X 

X 

X 

n. i 

N. 

fw 

N. 

c 

B 

B 


O 

B 

B 

B 

C 

a 

B 

B 

B 

B 

B 

s 

B 

B 

B 

B 

B 

B B 

B 

B 

BBS 

B 

B B B B 

B 

B 1 

B 

B 

B 


ORIGINAL PAGE 
COLOR PHOTOGRAPH 


WdSit iKmimui 


205 


y hi v 


PRECEDING PAGE BLANK NOT FILMED 



CONTOUR l EUELS 



CCOSOCOQCCSCCCCCCCI 

Sce = ccc = == = = cc=c: = ccccccc 
©xEuioinoifiomcme — — — — — 
dl N O N ui N O Nl,*. N c N x*. 

— — Px^MP^nrsrjnTTT 
doDcddddede 


ooooooocoocoo 
Eoopooqoooooo 

IT C X C X © X BlIClflBlflOl/lBinDlflOlflOineir 
rwcr^tn^ocsiinf^acsim^Ofsitn^ocsiin^-pcsiknN- 
, rxxxxxxxifiN*N-N-r^®®®®®Ciff>o»oooo 


e c s e e < 


C C C C I 


ICOCOOOO — — — — 


ORIGINAL PAGE 

COLOR PHOTOGRAPH 


207 


MENTION Al i r mi 


PRECEDING PAGE BLANK NOT FILMED 


Fig. 102 Lower Surface Pressure Contour, t-5.157. 
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Fig. 103 Lower Surface Pressure Contour, t-5.498. 
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Fig. 104 Upper Surface Pressure Contour, t-4.138. 
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Fig. 105 Upper Surface Pressure Contour, t-4.56. 
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Fig. 106 Upper Surface Pressure Contour, t-5.157. 
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Fig. 107 Upper Surface Pressure Contour, t-5.498. 
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